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We develop a quantum theory of electron confinement in metal nanofilms. The theory is used to 
compute the nonlinear response of the film to a static or low-frequency external electric field and to 
investigate the role of boundary conditions imposed on the metal surface. We find that the sign and 
magnitude of the nonlinear polarizability depends dramatically on the type of boundary condition 
used. 
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I. INTRODUCTION 



The recent explosion of growth in the field of nanoplas- 
monics^— have been a joint success of theory and exper- 
iment. Yet, in certain respects, theory is lagging be- 
hind. One profound theoretical question, which has not 
received adequate attention so far, is related to the ap- 
plicability of macroscopic theories. That is, the theory 
of plasmonic systems is almost exclusively based on the 
macroscopic Maxwell's equations, even though the sam- 
ples involved are, in some cases, only a few nanome- 
ters in size. The problem is compounded by the fact 
that plasmonic applications utilize highly conductive no- 
ble metals. In this case, the mean free path of the con- 
duction electrons, which can be significantly larger than 
interatomic distances, becomes of primary importance. 
Several different approaches to accounting for the finite- 
size and quantum effects in nanoparticles have been pro- 
posed^. However, finite-size-dependent electromagnetic 
nonlincarities have received relatively little attention. In 
this paper, we theoretically investigate finite-size effects 
in metallic nanofilms with an emphasis on electromag- 
netic nonlinearity and on the role of boundary conditions 
applied at the nanofilm surface. 

A qualitative understanding of finite-size effects in 
plasmonic systems can be gained from purely classical 
arguments. For example, it has been demonstrated that 
the electromagnetic response of a strongly coupled plas- 
monic system is dramatically altered by the effects of 
nonlocality when the smallest geometrical feature of the 
system (e.g., an interparticle gap) becomes comparable 
to a certain phcnomenological length scale, which char- 
acterizes the intrinsic nonlocality of metalZr— . A quali- 
tative theory of optical nonlinearity can also be derived 
from purely classical arguments^. 



However, a quantitatively accurate theory must be 
quantum mechanical. Unfortunately, a fully self- 
consistent and mathematically-tractable quantum model 
of p lasmonic systems is difficult to formulate. In Refs.lTTl. 
Il2i . a metal nanosphere was approximated by a degener- 
ate electron gas confined in a spherical, infinitely-high 
potential well. The conduction electrons were assumed 
to be driven by a time-harmonic and spatially-uniform 
electric field. This driving field is internal to the nanopar- 
ticle; relating it to the external (applied) field constitutes 
a separate problem-^. The model of Refs. [TTIIT2I has at- 
tracted considerable attention in the optics community. 
In particular, it predicts correctly the size dependence 
of the Drude relaxation constant. Although the model 
is mathematically tractable, it yields an expression for 
the third-order nonlinear polarizability, which involves an 
eight-fold series. This expression can be evaluated only 
approximately^. Recently, we have reduced analytically 
this expression (without resorting to any additional ap- 
proximations) to a five-fold series, which is amenable to 
direct numerical evaluation^. We then have used this 
result to study the size- and frequency dependence of the 
third-order nonlinear polarizability. 



Still, formulation of the model of Refs. [Tllfl^ in- 
volves two important assumptions, which arc difficult to 
avoid and whose effect is difficult to predict. First, the 
model assumes a spatially-uniform electric field inside the 
nanoparticle. In reality, this field is not spatially-uniform 
because the induced charge density is different from the 
infinitely-thin surface density of the Lorenz-Lorentz the- 
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ory, which was used in Refs. I12II13I and also by u: 
obtain a relation between the internal and the applied 
fields. A relation of this kind is essential for deriving 
any experimentally-measurable quantity. Unfortunately, 
a more fundamental approach to relating the internal and 
applied fields does not readily present itself, at least, not 
within the formalism of Refs. Illlfl2l . Second, the impo- 
sition of an infinite potential barrier at the nanoparticle 
surface has not been justified from first principles, even 
though this boundary condition is frequently employed 
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and can be defended by noting that a gas of nonintcract- 
ing electrons can not be stably bound by the Coulomb 
potential of a spatially-uniform, positively-charged back- 
ground. In other words, the potential barrier makes up 
for the neglect of the discrete nature of positively-charged 
ions. 

The two assumptions described above can be, in prin- 
ciple, avoided by using the density-functional theory 
(DFT). Within DFT, the exchange-correlation potential 
renders the system stable even if the ionic lattice is re- 
placed by jellium. Time dependent DFT (TDDFT) and 
the jellium model have been successfully used to com- 
pute the linear response of nanoparticle o 15 ' 16 . TDDFT 
and the random-phase approximation (RPA) have been 
used without resorting to jellium model (that is, with the 
full account of crystalline structure of metal) to study the 
effects of surface molecular adsorption on the dielectric 
losses in meta l 17 i 18 . Size- and shape-dependence of the 
imaginary part of the dielectric function of Ag nanopar- 
ticles has been also studied experimentall y 19 ! 20 . 

However, it is not clear whether the same DFT-based 
approaches are appropriate for obtaining nonlinear cor- 
rections to the polarizability. Indeed, the exchange- 
correlation potential is typically constructed for an in- 
finite and spatially-uniform system and is not expected 
to be accurate near boundaries, where the electron den- 
sity changes on small spatial scales. Yet, it can be argued 
that the nonlinear response of a nanoparticle is strongly 
influenced by the electron density near the boundary^ -. 
Further, if the jellium model is used, as is the case in this 
work, the binding potential of the positively-charged jel- 
lium (even with the account of exchange-correlation po- 
tential) is not as strong as that of discrete ionic lattice. 
In the latter case, the Coulomb potential approaches neg- 
ative infinity in the vicinity of the ion cores. Therefore, 
it is not evident whether the correct model should incor- 
porate an additional potential barrier at the surface of 
the nanoparticle to account phcnomcnologically for the 
reduced binding power of the jellium. Considering the 
above uncertainty, and the fact that the infinite poten- 
tial barrier of Ref. [l2lfl3l has been widely used in the 
optics literature, it seems desirable to investigate the in- 
fluence of the boundary conditions used on the obtained 
nonlinear corrections. 

To address the problem formulated above, we have 
performed DFT calculations for the simplest plasmonic 
system - a thin metallic film to which a static or 
low-frequency perpendicularly-polarized electric field is 
applied. We have investigated the influence of three 
different types of boundary conditions, including the 
rigid (infinite potential barrier) boundary condition, the 
Bardeen's boundary condition (in which the potential 
barrier is displaced from the metal surface) and free 
boundaries. We compute nonlinear corrections to the 
film polarizability by two different approaches - first, by 
the use of a perturbation theory, which is described in de- 
tail below and, second, by direct numerical application 
of the DFT (for control). For relatively strong applied 



fields, when the perturbation theory is not applicable, 
the nonlinear polarization is computed directly by DFT. 

The paper is organized as follows. In Sec. |TT] we de- 
scribe the physical model in more detail. In Sec. [TTTJ 
we write the basic equations of the DFT and specialize 
them to the thin film problem. The perturbation theory 
is described in Sec. IIVI Numerical results are reported 
in Sec. |V] Finally, Sec. IVII contains a summary and a 
discussion. 

II. THE PHYSICAL MODEL 

Throughout this paper we consider the following model 
system. A thin metallic film of width h is taken to occupy 
the spatial region — L/2 < x,y < L/2, —h/2 < z < h/2, 
where L h. The latter condition allows for the neglect 
of effects due to the edges of the slab. Accordingly, we as- 
sume that all physical quantities depend only on the vari- 
able z, so that the system is effectively one-dimensional. 
A static electric field of strength <g is applied to the slab 
in the z-direction. We emphasize that this is the external 
(applied) field; not the internal field of Refs. [Tllfl^ . 

We use the stabilized jellium model2i~— and the 
Gunnarsson-Lundqvist local-density approximation for 
the exchange-correlation potential 2 ^ (defined precisely 
below). At or near the film surface, we apply three differ- 
ent types of boundary conditions. In our first set of cal- 
culations, we apply the infinite potential wall condition of 
Refs. llflfl2l at the physical boundary of metal. We denote 
this type of boundary condition by "R" (in memory of 
S.G.RautianF^). In a second set of calculations, employ- 
ing what we refer to as the "B"-type boundary condition, 
we use the classical Bardeen model 2 ^— in which the po- 
tential wall is displaced from the physical boundary by 
the distance 



where k F is the Fermi wave number. The displacement 
Ab can be expressed in terms of the characteristic length 
scale of the problem, i, where £ 3 is the specific volume 
per conduction electron. Using the expressions Hk F = 
y/2m e E F and E F = {Z^f^K 2 /2m e P- , with E F being 
the Fermi energy and m e the electron mass, we find that 
A B = (9tt) 1 / 3 8- 1 ^ w 0.38^. Finally, in the third set 
of calculations, we do not use any additional potential 
barriers at or close to the metal surface. We denote this 
last boundary condition as type- "F" . 

In all these cases, we compute the electromagnetic re- 
sponse (polarization). The quantity of interest is the in- 
duced dipole moment per unit area of the slab 3 s , which 
is related to the charge density of the conduction elec- 
trons -pe(z), by 

9 = - J zp e (z)dz . (2) 

The quantity p e (z) > has been computed using stan- 
dard DFT theory using the three different types of 
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boundary conditions, which are described above. We 
have investigated the dependence of 3? on the strength of 
the applied field and on the film width. When nonlinear 
corrections to 2? are computed, the characteristic scale 
for the applied field strength is the atomic field 
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T- 



(3) 



The surface density of the dipolc moment can then be 
conveniently normalized by the quantity 



(4) 
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For a relatively small applied field, the induced polar- 
ization can be expanded in powers of S'/S' a ± according 
to 
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Note that the coefficients a.?., 04, etc., are identically zero 
in the above expansion due to the slab symmetry. 

The slab width can be parametrized by the number of 
atomic layers, M. Many metals of interest in plasmonics 
(including silver) have an fee lattice structure with four 
conduction electrons per unit cell. In this case, the slab 
width is h = Ma, where the lattice step a is related 
to I by a = 4 1 / 3 £. The slab then occupies the region 
-Ma/2 < z < Ma/2. In silver, I ps 0.26nm and a 
0.41nm. 



III. BASIC EQUATIONS 

The starting point in DFT is the eigenproble m 28 i 29 : 

h 2 



2m, 



■V J + U eS (z) 



Vv(r) = EMv) , (6) 



where the index p labels energy eigenstates, m e is the 
electron mass and 



U eS (z) = U K (z) + U KC (z) + V(z) 



(7) 



is the effective potential consisting of the Hartree term 
Un(z), the exchange-correlation potential C/ xc (z), and of 
the interaction potential 



V(z) 



-ezS 



(8) 



which describes interaction of the electrons with the ap- 
plied field. 

We adopt the stabilized jellium model 2 ^— , according 
to which the positively-charged ions form a spatially- 
uniform charge density Pi(z), such that pi{z) = e/£ 3 if 
\z\ < h/2 and Pi(z) = otherwise. In the stabilized jel- 
lium model, a constant potential (in our case, negative- 
valued) is added to U e g(z) inside the spatial region oc- 
cupied by jellium. This constant potential is chosen so 



as to make the metal mechanically stable at its observed 
valence electron density. As we found, the difference in 
the results obtained using the usual and the stabilized 
jellium models is insignificant (it is absent for the It- 
type boundary condition). On the other hand, the use of 
the stabilized jellium model removes certain anomalies of 
the standard jellium model, such as negative surface en- 
ergy 2 ^. Also, in the case of a silver film, we have used the 
stabilized jellium to compute the work function, which 
we find to be 3.8eV for the thickest film considered (with 
M = 32 atomic layers). This result can be compared 
to the experimental measurements in silver— [4.2eV for 
the (100) face]. We have obtained a somewhat less ac- 
curate prediction without jellium stabilization (3.5eV). 
Therefore, the stabilized jellium model is used in all cal- 
culations reported below. 

The Hartree potential is given by 



Un{z) = 2ne [ Pi (z') - p e {z') 



>\dz' , (9) 



where —p e {z) is the density of charge associated with the 
conduction electrons. Both densities are normalized by 
the condition J pi{z)dz = J p e (z)dz = eh/£ 3 . 

The Gunnarsson-Lundqvist local-density approxima- 
tion for the exchange-correlation potential, which we use 
in this work, is of the following form: 



U xc (z) 



Cx 



R 



C c — lull 

(IB 



A 



tip, 
R 



(10) 



Here R = {Se/Airp,,) 1 ^, a B = h 2 /m e e 2 is the Bohr ra- 
dius, and the dimensionlcss constants are C x = 0.61, 
C c = 0.033, and A = 11.4. 

Assuming periodic boundary condition at x = ±L/2 
and y = ±L/2, we can write the eigenfunctions of ([6]) in 
the form 



exp 



1 — (lx 
L 



my) 



(11) 



where <f> n (z) satisfies the equation 

h 2 



2m, 



In Eq. (|11|) . we have replaced the index \i by the triplet 
of quantum numbers (l,m,n). In what follows, we view 
p as a composite index. For example, summation over p 
entails summation over the three quantum numbers I, m 
and n. The energy eigenvalues are given by 



Eimn = (27rfr) : 



I 2 + m 2 
2m P L 2 



(13) 



The ground-state electron density is given by the expres- 
sion 

Pe (z) = 2 e ]Te(£ F -i^)l^(r)| 2 

Tip 

= eY / W n \Mz)\ 2 ■ (14) 
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Here the factor of 2 in the first equality in (1141) accounts 
for the electron spin, Ep is the Fermi energy, np is the 
maximum quantum number n for which there exist en- 
ergy levels Ei mn below the Fermi surface, and W n is a 
statistical weight given by 



m e (E F - e n ) 
irh 2 



(15) 



Note that the Fermi energy E-p for the model consid- 
ered in this paper is close to the macroscopic limit of 
(h 2 /2m e )(3ir 2 /£ 3 ) 2 / 3 . In all calculations shown in the 
paper, we have computed E-p and the related quantity 
np by ordering all energy levels E\ mn . It should be also 
emphasized that the quantity np is determined with the 
account of all degeneracies, including those due to the 
spin, so that the second equality in ([M]) does not contain 
the factor of two. 

Once p e {z) is found, the dipole moment per unit area 
of the slab 8P can be computed according to @. 



IV. PERTURBATION THEORY 

A perturbative solution to the eigenproblem (|T2"j) is 
obtained by expanding the quantities 4> n {z), e n , W n , 
U e g(z), p e (z) into power series in the variable x = S I S^- 
For example, we write 



(16) 



Upon substitution of these expansions in the original 
eigenproblem ([T2|) . we obtain the following relations (to 
third order in x): 



s = 0: 



s = 1 



(i) -u 77PLU0) = Jo)Ai) , ti)A0) 



s = 2 



(17a) 
(17b) 



JX 



. ( 17c ) 



s = 3 : 



e W^)+ e f0« +e (3)^(o). (i 7d ) 



Here 



H = -J_ 5 2 + C/ (0) + t7 (0) 



(18) 



2m e ~ 2 1 ff xc 
is the unperturbed Hamiltonian and 

/oo 
-00 

+ ^(^W(,) + yW(z), (19) 



where 



V {0 \z) =0 , 
0%) 

^ 2 >(z) 



(20a) 
(20b) 

(20c) 



V^(z) = U^(z)p^{z)^ 2 \z) + \u™\pV(z)]*. (20d) 

The operators U^ c , U" c , and C/" c ' are the first, second, and 
third functional derivatives of the exchange-correlation 
potential /7 XC with respect to p e evaluated at p e = 

In deriving (fH?)) . ([2"0")l , we have taken into account the 
implicit dependence of U e s on the expansion parameter 
x, which stems from the dependence of p e on x. The first 
four coefficients in the expansion of p e can be obtained 
using equations (fT4]) and p5[) . which results in 



2W. 



(0)^,(0)^(1) 



n T'n rn 



=e 



,(3) 



E (^ 0) ) 

n=l - 

"F r 2 

E ^i 2) (^ 0) ) +2^ 0) 4 0) 4 5 

n=l L 

(<#>) 

; eK o) (4 o) 4 3) +4 i) 4 2) )+ 



(21a) 



(21b) 



2VFi 2 ) 



n rn rn 



Here 



wi s) 



E< 

,k=i 



ripe 



(21c) 



(22) 



are the coefficients in the expansion of the statistical 
weights and 



w 



imp 



h 2 



(23) 



We note that the integer number np, which was defined 
(after Eq. 1141) as the maximum quantum number n for 
which there exist energy levels Ei mn below the Fermi 
surface, is independent of the applied field § for S < 
10^ at . Therefore, in the perturbation theory developed 
here, it is sufficient to compute np at zero applied field. 

The procedure of constructing the perturbation series 
for the dipole moment density, starts from solving 
(|17a[) numerically. The eigenvalues el are ordered so 
that e n < e„+i for n = 1, 2, .... It follows from symmetry 
considerations that 



J(0) 
JO) 



i-z) = {-ir^^\z) 

(0)/ 



P^(-z) = p?>(z) 



(24a) 
(24b) 
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These relations have been strictly enforced in the numer- 
ical procedures we have implemented. 

We seek the solutions to (|17l) in the following form: 



(*)• 



(25) 



Here the unperturbed basis functions \4> k ) (and the en- 
ergy levels €n^) are assumed to be known; they are deter- 
mined by solving (|17a[) numerically. In the zeroth order 

(s = 0), we have trivially C^. = 5 n k- In higher expan- 
sion orders, we determine by substituting (f25|) into 
(fT7| . For s 
equations: 



1,2,3, this results in the following three 



«') E i 



JO) 



(2) 



> = E 



i^ 0) > 



.(0) JO) 



(26a) 



(26b) 



(0)\ 



.(0) 



JO) 



= E 

+ (4°Ve ( ff ) i^ ) > + (4 0) i^ ) i^ ) )-(4 0) i^ 1) ) 

x f^l^l^ + ^lt^l^})! • (26c) 



Regarding the set of equations (f2l))) , we note the follow- 
ing. First, the vectors ) must be computed recur- 
sively starting with s = 1. Thus, we compute \4>n^} ac- 
cording to (|26a|) . The result is substituted into (|26b[) . 

which allows one to compute \4> n ), and so forth. How- 
ever, unlike in the ordinary perturbation theory, the 
operators in the right-hand sides of (|2"6")l arc still 
unknown and must be determined separately. Indeed, 
U^fi(z) depends on the functions p^ s \z) according to 
(fH?|) and the latter depend on the functions 4>n ' (z) with 



s = 0,1, ... ,s according to (|2TJ) . and some of these func- 
tions cf>n \z) have not yet been computed, even at s = 1. 
We therefore must consider equations (fTi?|) . (|2"Tf and (|2"o) 
as a coupled set of equations, which must be solved self- 
consistently. 

To compute the matrix elements appearing in the 
right-hand sides of (|2"6"|) . we use the following procedure. 
We start with s = 1 and define a (yet unknown) quantity 



A 



(i) - 



(i) 



In particular, we have e„ 



(27) 



w n 1] , we can express pe (z) in terms as follows: 



n=l 



U=l 



WE 



(o) _ (o) 



(4 0) W) 2 + 

(28) 



4%) 



We then substitute (|28|) into (fT9|) (in which we must spe- 
cialize to the case s = 1) and compute the expectation 
value of U^g between the unperturbed states (4>n^ | and 



k=l l^k 



to obtain the following set of linear equations: 

fc=i \i=i / 
x(4 0) l4° fc ) l€ ) >=< 1 2- (29) 



This set must be solved with respect to the unknown 
quantities A^m ■ The matrix L is defined by the following 
relations, which contain only the known quantities: 



2W, 



(0) \ <Pn |G fc; 1 0m ) 
JO) JO) 



(30) 



G<°>(*)=27re a / \z-^(04 0) (0^ 



K c {z)^\z)<t>?>{z 



(0), 



(31) 



and the right-hand side of the equation is determined by 
the following expression: 



(32) 



The computational complexity of solving (|29[) can be 
reduced by making the following observation. It fol- 



(i) 

m+2q,m 



= o for all 



lows from (|20b| and (J22) that R, 
q = 0,±1,±2, ... and m,m + 2q > 1. Using $fl$ and 
the relation 



(33) 



we find that the set of equations (|29|) splits into two un- 
coupled subsystems. The first subsystem is homogeneous 



and contains the matrix elements of the form A 



(i) 



m-t-2g,m ' 

The only solution to this subsystem is trivial. Thus, we 



have A 



(i) 

m+2q,m 



0, which also implies that the last term 



in the left-hand side of ([29| vanishes and eh^ = Ann = 0. 



The other subsystem is of the form 



= R 



(i) 



A W , rfc + 2p+l,fc A (l) 

' LA -ro+2q+l,m ~r / j / j J - / m+2q+l,m'- y k+2v+l.k ~ -"*m+2g+l,m 
k=l p 



(34) 

By substituting the with a nonzero right-hand side. Eq. (jMJ) must be solved 

(i) 



expansion (|26aj) for \(jr n ) into (|21a[) and using (|2"2"j) for numerically with respect to the unknowns A m+2q+1 
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In addition to the simplification outlined above, the fol- 
lowing parity relations hold: 



<t><t\-z) = (-ir<i>W{z), 

pP{-z) = -p ( i-\z), 



(35a) 
(35b) 
(35c) 



After solving (|34|) numerically, one can find \<pn^} from 
([26a]) . and also p ( e\z) and U^(z) from (|2Ta|) and (fT9|) . 
respectively. 

Next, we take s = 2 and start again with computing 

(2) 

the matrix elements A kn , which are defined in a manner 



similar to (|2"7|) : 



(36) 



Once are found, \4>n^) can be obtained using (|26b>[) 
while the second-order correction to the energy of the 
n-th level is given by 



(37) 



Note that the terms \U^ and (^l^) in 

(|26b[) are first-order quantities, which have already been 
computed. 

After substituting (|26b|) and (J22J) into (|21b|) . we ar- 
rive at the expression for pe (z), which is similar to (j2"8")l 
and is omitted here. The only unknown quantities in 

(2) 

this expression are the matrix elements Aj.„. We then 
substitute this expression into (|19p , in which we special- 
ize to the case s = 2, and take the expectation between 
the unperturbed states. This yields the following set of 
equations: 



k=l l^k 

x (</>!°)|G (0) u(°> 



'kk \Ym 




(2) 



(38) 



We see that (l3l) has the 



Here the unknowns are A 
same matrix as (|29[) but a different right-hand side, viz, 

"F 

rip 

-A*<^|Gffi|^>-<^|GS|^>l 



2e z 



(39) 



where 



G^(*) = 27re 2 



UUzM\z)] 2 



and 



A, = 



E< 

i=i 



'i 



- n F (0 fc |U cff |0 fe 



off I s^z 



(41) 



? (2) 



All quantities, which enter the definition of R nm , are 
determined by in the first-order. As follows from (|24al) . 

(22), and dSSJ), fl^ 2 g+l,m = for a11 3 = °> ±1> ±2 > • • • 
and m, m + 2g + 1 > 1. Thus, unlike in the case s = 1, 

(2) 

we now have a homogeneous subsystem for A^i^g+i m , 
which has only the trivial solution, while the quantities 

(2) 

A y m ' +2qm satisfy the following inhomogeneous subsystem: 



A 



(2) 



m+2q,m ' J ^m+2q,m'- :s 'l+2p,l 



(2) 



1=1 p^o 



lip / np 



I? 



fc=i \;=i 

(2) 



(42) 



(21 

We then solve (|42|) with respect to A^^_ 2(J m numerically 
and use the result to compute the second-order correc- 

(2) (2) (2) (2) 

tions, that is, the quantities \4> n ), pi , , Wn , and 
el 2) by using equations ([26b]) . Iplb]) . (15]). and ([3T|). 
respectively. At the second order, the following parity 
relations hold: 



4 2 )(-z) = (-ir-vi 2) w, 



pi 2) (-^) = p™ 



(2), 



(2), 



(43a) 
(43b) 
(43c) 



Finally, in the case s = 3, all quantities in the right- 
hand side of (|26c|) arc known except for 



(44) 



Again, once A^ are found, \4>^) can be obtained from 
(|26cp while the third-order corrections to the energies are 
given by 

=A(f) + mU^) - eW(0\^) 

+ (0W|C/ c (2 Vi 1) > • (45) 



It follows from the parity relations (|24a[) . (|35l) . and (|43|) 
that all terms in the right-hand side of (|4"5|) vanish ex- 
cept for A^ 3 2, so that we have e„ = A„„. In fact, we 

(3) (3) 

will see below that e„ = A„„ = 0. Similarly to the pro- 
cedure outlined above for the s = 1 and s = 2 cases, we 
substitute (j) ( n ] from ([2"6"c]l and from $12) into ([2"Tc]) . 
obtain an expression for pe (-2)7 in which the only 1.111- 

(3) 

known quantities are , and substitute the result into 

(40) 

(fT9")l . where we specialize now to the case s = 3. Finally, 
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we compute the same expectations as before and obtain 
the set of equations 



Tip / Tip 



a< 3 > + E E J&Ai? + - E E < 3) - »f A i 

fe=i fe=i \/=i 

X(^)| G (0)|^(0) )=J R(3) 1 



(46) 



(3) 

with respect to the unknowns A„m- The matrix is the 
same as before, while the right-hand side is given by the 
following relations: 



R (3) 

- 1 ' II ill 



>Y, W > 

k=l 



(0) 



E^1 0) r»> 



k=l 

+ l(40)| O (l) (2) | ^) ) 

-^mU^?m- (47a) 



^ c (z)0i 1} (z)0i O) (z) , 



(48) 



G$(z) = 27re 2 / Iz-d^^O^m 

*/ — OO 



(49) 



where 



JO) (0) 



(3) 

Similarly to the case s = 1, we have R m +2q m = f° r a ^ 
(7 = 0, ±1, ±2, . . . and m,m + 2q > 1. This follows from 
the parity relations (|24a[) . ([33]). ([35]). and (|i3"]). Conse- 



(50) 



(3) 

m+2g,m 



and el 3) = A^ = 0. The 
,0) 



qucntly, we have A 

matrix elements of the form A]^' +2q+1 m are determined 
from the following subset: 

Tip 

A (3) jk+2p+l,k a (3) _ R (3) 

AA m+2o+l,m~ 1 "/ , / y ■ L, m+2g+l.m^ A fc+2p+l.fc — ' n -m+2g+l,m ■ 
fe=l p 

(51) 

As in the s = 1 case, the following parity relations hold: 



pf\-z) = -pf\z) , 
ujg(-z) = -Ug\z). 



(52a) 
(52b) 
(52c) 



By solving (|5T|) we can find, in particular, the expansion 

coefficient p^\z) for the negative charge density from 
(f2"Tc|) using (|26"c|) and ([22]) at s = 3. The expansion 
coefficients a\ and 03 in (|5|) arc then obtained from 



47T 



^ s) (z)dz . 



(53) 



at J —00 



It can be seen from ()43b[) that 0L2 = 0. 



V. RESULTS 

We have performed computations for varying film 
width and different boundary conditions. We have found 
that, in the case of R- and B-typc boundary conditions, 
the effect of including the exchange-correlation poten- 
tial is relatively minor. However, in the case of the F- 
type boundary condition, the exchange-correlation po- 
tential must be included in order to stabilize the con- 
duction electrons. For the purposes of a fair comparison, 
we show the results of R- and B-type simulations with- 
out the exchange-correlation potential, except in Fig. 3 
below, where the results with and without the exchange- 
correlation potential are compared. It will be shown that 
the R- and B-type boundary conditions produce a nega- 
tive nonlinear correction to the polarizability and satura- 
tion effects, in agreement with Refs. [l0l - [l2l . However, the 
F-type boundary condition results in a positive nonlin- 
ear correction whose magnitude is a few hundred times 
larger than that in the case of R- or B-type boundary 
conditions. In all cases, the emergence of macroscopic 
(bulk) behavior becomes evident in relatively wide films. 

The eigenproblem (|12[) was solved algebraically by 
discretizing the differential equation in the interval 
-z m ax/2 < z < z m ax/2. Here z max = h+12a = (M+12)a 
for the F-type boundary condition, z max = h + 2A^ for 
the B-type boundary condition (As is defined in (p}), 
and z max = h for the R-type boundary condition. Recall 
that the jellium (the physical slab) is contained in the 
region — h/2 < z < h/2, where h = Ma. Central differ- 
ences with 20 discrete points per the lattice unit a have 
been used and convergence was verified by doubling this 
number. 

In Fig. 1, we plot as a function of the applied field, 
for different widths of the slab and for different types 
of boundary conditions. In the case of the R- or B-type 
boundary conditions, the system is stable independent of 
the strength of the applied field. In the case of F-type 
boundary condition, we have observed an instability of 
the charge density for $ > 0.1<oat- This instability can 
be explained by the effect of tunneling, which can result 
in significant charge accumulation in a biased potential 
over long periods of time, or after many DFT iterations. 
The data points affected by this instability are not dis- 
played in Fig. 1. It can be seen that the F-type boundary 
condition tends to increase & (compared to the macro- 
scopic limit) while the B- or R-typc boundary conditions 
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FIG. 1: The dipole moment density per unit area, SP, as a 
function of the applied field, £, for different types of bound- 
ary conditions and for films consisting of different numbers 
of atomic layers, M. The continuous line corresponds to the 
"bulk" (macroscopic) result and the centered symbols repre- 
sent the results obtained from the DFT. In panel (a), three dif- 
ferent values of M are used, as labeled. All data points above 
the continuous line correspond to F-type boundary conditions 
and the points below the line correspond to R-type boundary 
conditions; B-type boundary conditions are not used in this 
panel. In panel (b), only M = 2 is used. The two sets of 
centered symbols below the continuous line in (b) correspond 
to R- and B-type boundary conditions, as labeled. In both 
panels, the data points for F-type boundary conditions (above 
the continuous line) terminate at S /<f a t = 0.87; computations 
for larger values of £ (with F-type boundary conditions) are 
affected by the numerical instability discussed in the text. 



tend to decrease 3P . For all types of boundary condi- 
tions, deviations from the macroscopic result are signifi- 
cant when M = 2 and M = 8 but small when M = 32. 
In the case M = 2 (Fig. lb), the B-type curve lies above 
the R-typc curve; a similar result was obtained for other 
values of M (data not shown). Physically, the behav- 
ior illustrated in Fig. 1 can be understood as the result 



of electron spillove r 31 ' 32 (in the case of F-type boundary 
condition) or as the combined action of the finite charge 
density of the jellium and of the uncertainty principle (in 
the case of B- or R-type boundary conditions). 

Although the deviation of from the macroscopic 
result is obvious in Fig. 1, all curves shown in this figure 
appear to be linear. To visualize the deviations from 
linearity, we have computed the nonlinear contribution 
to the dipole moment density according to 



nonl 



ha.\ „ 
& -§ 

47T 



(54) 



The result is plotted in Fig. 2. It can be seen that the 
correction is negative for R- and B-type boundary condi- 
tions. For the F-type boundary condition, the correction 
is positive and about 200 times larger in magnitude. An 
interesting effect can be seen in Fig. 2b. Namely, the B- 
type boundary condition produces a nonlinear correction 
of a larger magnitude compared to the R-type boundary 
condition. This is somewhat unexpected since the data 
of Fig. 1 suggest that the B-type curve is closer to the 
macroscopic asymptote. Moreover, we have discovered 
a non-monotonic dependence of the expansion coefficient 
a 3 in ([5]) on the displacement parameter A, as is dis- 
cussed below. 

Recall that the B-type boundary condition involves a 
displacement of the rigid potential wall from the surface 
of the metal by the distance A = Ab ~ 0.381 We 
can, however, view A as a free parameter. In the case 
A = 0, we recover the R-type boundary condition, in 
the case A = oo - the F-type boundary condition, and 
A = Ab corresponds to Bardeen's model. In Fig. 3, we 
plot the expansion coefficient as a function of A for 
M = 2. In this figure, we show the data obtained both 
with and without the exchange-correlation potential. We 
find that the surprising non-monotonic dependence is ob- 
served in both cases. For larger values of A, the red curve 
in Fig. 3 (with exchange-correlation potential included) 
rapidly grows and saturates at the level of CV3 « 0.1 for 
A > 8Ab (data not shown). The latter result exactly 
corresponds to the one obtained with the F-type bound- 
ary condition. Note that ol\ is almost independent of 
A. Qualitatively the same results have been obtained for 
M = 8 and M = 32. 

Next, in Fig. 4, we plot for the B-type boundary 
condition (with A = Ab) in a very large interval of £*, 
up to §j§ax = 10 3 . Of course, an applied static elec- 
tric field of this magnitude is not achievable in practice. 
However, the situation can be more experimentally fa- 
vorable in the case of quasistatic fields. Although we 
do not consider this case directly, it is known that the 
internal field enhancement factor due to plasmon reso- 
nances is of the order of uip/j, where uj v is the plasma 
frequency and 7 is the Drude relaxation constant. This 
factor can be as large as ~ 500 in the case of silver, and 
it enters the nonlinear correction to the polarizability in 
the fourth powe r 13 ' 14 . Note that the quasistatic approx- 
imation (known in the context of DFT as the adiabatic 
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FIG. 4: # as a function of the applied field, S, for B-type 
boundary condition and for different numbers of atomic lay- 
ers, M. Centered symbols correspond to DFT results and 
continuous lines to the analytical expression (|55|) . 



approximation) is applicable as long as hcj/c <C 1, which 
easily holds even in the visible spectral range. Also, the 
electric field intensity in very short laser pulses can be 
of the order of or higher than the atomic field, and the 
related physics has attracted considerable recent atten- 
tion^. 

In Fig. 4, we also compare the DFT calculations with 
the expression 



FIG. 2: J^noni as a function of the applied field, S . The 
same convention for encoding the different types of boundary 
conditions and the different numbers of atomic layers as in 
Fig. 1 is used. Additionally, different numerical normalization 
factors JY have been used for different types of boundary 
conditions, as labeled. 
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FIG. 3: Dependence of a 3 on A for M — 2, with and without 
accounting for the exchange-correlation potential (ECP). 
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(55) 



which was derived in Ref. [l0| using purely classical ar- 
guments. It can be seen that (|55|) is surprisingly accu- 
rate for M = 8 and especially for M = 32. This may 
seem unexpected because (|55p contains a nonanalytic- 
ity of the form &\$\, while the expansion ([5]) represents 
a real analytic function. This discrepancy is resolved by 
noting that ([5]) has a finite radius of convergence and that 
an expansion of this type can not, in principle, capture 
the saturation phenomena illustrated in Fig. 3. On the 
other hand, numerical DFT calculations can be carried 
out whether or not (O converges. 

Finally, we investigate the dependence of the coeffi- 
cients ai, as on the number of atomic layers M for 
all three types of boundary conditions. The results are 
shown in Fig. 5. It can be seen that, in all cases, the de- 
pendence is monotonic. Comparing the results for R- and 
B-type boundary conditions, we reconfirm the trend that 
has been already noted, namely, that B-type boundary 
conditions produce a smaller finite-size correction to ai 
(compared to R-typc boundary conditions) but a larger 
nonlinear response. 
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FIG. 5: Coefficients ai (a) and Q3 (b) as functions of the num- 
ber of atomic layers, M, for three different types of bound- 
ary conditions. Different numerical normalization factors jV 
have been used for different types of boundary conditions, as 
labeled. 



VI. SUMMARY AND DISCUSSION 

We have studied theoretically and numerically polar- 
ization of a thin silver film under perpendicularly ap- 



plied low-frequency external electric field. Three differ- 
ent boundary conditions have been applied at the film 
surface. It was shown that the sign and magnitude of 
the nonlinear correction to the film polarizability de- 
pends dramatically on the type of boundary condition 
used. Since all theories involved contain approximations, 
only comparison with experiment can determine which 
boundary condition is physically correct. 

An obvious shortcoming of the calculations reported 
herein is that they are carried out for static fields. How- 
ever, the results can be extended to finite frequencies, 
as long as relaxation and resonance phenomena are not 
taken into consideration, that is, if the frequency is far 
below the lowest plasmon resonance of the system. In 
practice, this means that the theory can be applied up to 
THz frequencies. Thus our results are amenable to exper- 
imental verification. A possible experimental test could 
be a measurement of the sign of the real part of the non- 
linear susceptibility f° r a suspension of silver nan- 
odisks at the excitation frequency ~ 10GHz. The disk 
thickness should be much smaller than the skin depth, 
S w 0.2/iin in this example. Previous experimental mea- 
surements of x*- 3 -* were largely confined to the optical and 
near-IR spectral regions 3 ^—, where the sign of Re[x (3) ] 
can depend on frequency due to the effects of plasmon 
resonances. 

It seems possible to further extend our theory to op- 
tical frequencies by utilizing the quasistatic approxima- 
tion, which is known as the adiabatic approximation in 
the context of DFT—. In this approximation, all po- 
tentials are computed using instantaneous values of the 
density, for example, by writing for the Hartree interac- 
tion potential Un[p e ](r,t) = Un[p e (r,t)}, and similarly 
for other functionals. This corresponds to neglecting the 
effects of retardation in the electromagnetic interaction 
and is adequate as long as huj/c <C 1. Note that plas- 
mon resonances and relaxation phenomena can be taken 
into consideration within quasistatics. However, time- 
dependent DFT is still relatively unexplored, although 
some promising results have been obtained 3 -^. 
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